Nonlinear optical properties of Te02 crystalline phases from first principles 
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We have computed second and third nonlinear optical susceptibilities of two crystalline bulk 
tellurium oxide polymorphs: Q-Te02 (the most stable crystalline bulk phase) and 7-Te02 (the 
crystalline phase that ressembles the more to the glass phase) . Third order nonlinear susceptibilities 
of the crystalline phases are two orders of magnitude larger than Q-Si02 cristoballite, thus extending 
the experimental observations on glasses to the case of crystalline compounds. While the electronic 
lone pairs of Te contribute to those large values, a full explanation of the anisotropy of the third 
order susceptibility tensor requires a detailed analysis of the structure, in particular the presence 
of helical chains, that seems to be linked to cooperative non-local polarizabilty effects. Our results 
demonstrate that first-principles simulations are a powerful predictive tool to estimate nonlinear 
optical susceptibilitites of materials. 

PACS numbers: 42.65.An, 78.20.Jq, 61.50.Ah 
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I. INTRODUCTION 

Tellurium oxide glasses arouse lots of interest in the 
field of nonlinear opties (NLO) sinee their unusual non- 
linear optical indices have been noticed. The third or- 
der optical susceptibilities [x^'^^] exhibited by pure Te02 
glasses, of the order of 14 x 10"^"^ esu^ are indeed 
among the highest observed for oxide glasses (50 times 
larger than in pure silica glasses) and they thus are of 
great interest in both fundamental science and techno- 
logical applications as optical modulators and frequency 
converters ji"— 

The origin of these high values is not fully established 
yet. Using a combination of experimental techniques (in- 
terferometric measurements) and ab initio calculations 
(within the restricted Hartree-Fock scheme), some au- 
thors^ suggested that the highly polarizable 5s^ elec- 
tronic lone pair of Tef^^^^ could be responsible for the high 
NLO indices. Other theoretical works;^ based on the 
use of hybrid functionals within the density functional 
theory (DFT), reinforced this idea by demonstrating the 
importance of the Te'^^^ lone pair on the hyperpolariz- 
abilities of isolated Te04 and TeOs structural units. This 
conclusion was extrapolated, by extension, to the case of 
Te02 based glasses. However, another series of theo- 
retical studies )2rJ^ also carried out in the framework of 
DFT with hybrid functionals on (X02)n (X = Si or Te) 
polymer clusters of different shapes (chains, rings and 
cage geometries) and sizes (monitored by the number n 
of XO2 units) suggested another origin for the unusually 
high values of the NLO susceptibilities, highlighting the 
relevance of the structural features themselves, in par- 
ticular how the structural blocks (i. e. Te02 units) are 
linked together. Only one type of such molecules, the lin- 
ear chains, seem to be capable of realistically reproduc- 
ing the high hypersusceptibility values for the tellurium 
oxides. This was attributed to an exceptionally strong 
nonlocality of the electronic polarization in these chains. 



much more important for the Te than for the Si oxides. 

Nevertheless, previous works were based on hypotheti- 
cal fragments that were supposed to be likely found in the 
glass. The question about the high nonlinear susceptibil- 
ity in the solid phases was not addressed. Clearly, further 
studies are needed to achieve a deeper understanding in 
the origin of these properties and notably to gauge the 
relative importance of the lone pair versus the structural 
features. A new way to tackle this problem is to treat the 
case of Te02-based crystalline compounds. This would 
prevent the recourse to hypothetical structural fragments 
to feed the first-principles calculations. Besides, it would 
allow to study how the anisotropic nature in the crys- 
talline phases translates into the variations of dielectric 
susceptibilities with crystalline directions. 

Unfortunately, the situation for crystals is different 
than for molecules and two main problems arise in the 
first-principles calculation of the hypersusceptibilities. 
The optical susceptibilities are derivatives of the bulk 
macroscopic polarization with respect to electric field. If 
the polarization can easily be expressed in terms of the 
charge distribution for molecules (finite systems), it can 
not be obtained that way for crystals (infinite systems 
treated periodically). The polarization in a periodic sys- 
tem would indeed depend on the choice of the unit celli^i 
Solutions arised in the early 1990s and are often referred 
to as the "modern theory of polarization" The basic 
idea is to consider the change in polarization^ of a crys- 
tal as it undergoes some slow change, e.g. a slow dis- 
placement of one sublattice relative to the others, and 
relate it to the current that flows during this adiabatic 
evolution of the systemJ^ The second problem lies in the 
nature of the applied electric field that is macroscopic. 
The scalar potential of a macroscopic homogeneous elec- 
tric field is non-periodic (so the Bloch theorem does not 
apply) , and unbounded from below (so the energy of the 
system can always be lowered transferring electrons to re- 
gions sufficiently far away, hampering the aplicability of 
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traditional variational methods) fl^ The first approach to 
circumvent this problem in first-principles simulations, 
due to Kunc and Resta,— was to consider "sawtooth" 
potentials in a supercell. We have previously tested this 
scheme, as implemented in the CrystalOG program^ii 
on the computation of hypersusceptibilities of Te02 crys- 
talline oxides.— Although this method gives satisfactory 
results, it requires defining a supercell for keeping the 
periodicity along the applied field direction. The dimen- 
sion of the studied system is thus very quickly limiting 
in terms of expansiveness in time and computational re- 
quirements. A more recent variational alternative, firmly 
rooted on the modern theory of polarization, was due to 
Souza, Ihiguez and Vanderbilt.— It is based on the 
minimization of a electric enthalpy functional with re- 
spect to a set of polarized Bloch functions, thus includ- 
ing the effect of the electric field directly inside the unit 
cell. This approach, recently implemented in the Siesta 
method^Si^ is the one used in the present work. 

In this paper we compute the second and third or- 
der optical susceptibility tensors in two bulk Te02 poly- 
morphs: the a-Te02 phase (known as paratellurite) that 
is the most stable one, and the 7-Te02 phase, whose 
structure ressembles the more to the glass. The estima- 
tions of the nonlinear susceptibility data provided in the 
present work intend to fill the gap in the reported values 
of these quantities. Unfortunately, there is a cruel lack 
of experimental nonlinear susceptibility data for crys- 
talline phases, mainly due to the difficulty of growing suf- 
ficiently large single crystals. To our knowledge, only the 
X^^^ susceptibility tensor elements for the a-Te02 phase 
has been measured, while there are no experimental val- 
ues of the x^'^'' tensor elements for any crystalline phase. 
In addition, third order susceptibility tensor of a-Si02- 
cristobalite, which is structurally similar to q;-Tc02, is 
computed for comparison. 

The rest of the paper is organized as follows. After 
presenting the computational details in Sec. HIl and the 
structure characteristics of the different polymorphs in 
Sec. mil we describe the methodology used to compute 
the nonlinear susceptibilities in Sec. IIVI Second-order 
susceptibility values are then calculated in Sec. IV Al and 
compared to experimental results in order to test the va- 
lidity and the limitations of the method. Finally, the 
method is used as a predictive tool through the calcula- 
tion of the third order susceptibilities in Sec. IV Bl and 
clues are given for exploring relevant features responsible 
for large variations of the dielectric susceptibilities with 
crystalline directions. 



II. COMPUTATIONAL DETAILS 

We have carried out density functional first-principles 
simulations based on a numerical atomic orbital method 
as implemented in the Siesta code.— All the calcu- 
lations have been carried out within the Generalized 
Gradient Approximation (GGA), using the functional 



parametrized by Perdew. Burke and Ernzerhof (PBE)^ 
to simulate the electronic exchange and correlation. 

Core electrons were replaced by ab initio norm con- 
serving pseudopotentials, generated using the Troullier- 
Martins scheme^^ in the Kleinman-Bylander fully non- 
local separable representatioujiSi The 5s and bp electrons 
of Te, 2s and 2p electrons of O, and 3s and 3p elec- 
trons of Si were considered as valence electrons and ex- 
plicitly included in the simulations. In order to avoid 
the spiky oscillations close to the nucleus that oftenly 
appear in GGA-gcnerated pseudopotentials, we have in- 
cluded small partial core corrections^^ for all the atoms. 
Te pseudopotential was generated scalar relativistically. 
The reference configuration, cutoff radii for each angular 
momentum shell, and the matching radius between the 
full core charge density and the partial core charge den- 
sity for the non-linear-core-corrections (NLCC) for the 
pseudopotentials used in the present work can be found 
in Table m 

The one-electron Kohn Sham eigenvalues were ex- 
panded in a basis of strictly localized^ numerical atomic 
orbitals322 The size of the basis set chosen was double- 
C plus polarization for the valence states of all the atoms. 
All the parameters that define the shape and the range 
of the basis functions were obtained by a variational op- 
timization of the energy in the a-cristobalite polymorph 
of Si02, and of the enthalpy (with a pressure P = 0.2 
GPa) in the a-phase of Te02, following the recipes given 
in Refs. [2^ and [2^. In both cases the optimization of 
the basis set was performed at the experimental lattice 
parameters and internal positions taken from Rcf. ,30, for 
Q!-Te02 and from Ref. |3l| for a-cristobalite. 

The electronic density, Hartree, and exchange correla- 
tion potentials, as well as the corresponding matrix el- 
ements between the basis orbitals, were calculated in a 
uniform real space grid. An equivalent plane wave cut- 
off of 400 Ry was used to represent the charge density. 
During the geometry optimizations, we used a 6 x 6 x 6 
Monkhorst-Pack mesh^ for all the Brillouin zone integra- 
tions. The macroscopic polarization and its derivatives 
with respect to an external electric field, depend highly 
on the number of fc-points used. To quantify this depen- 
dence, we have refined the Monkhorst-Pack meshes and 
followed the evolution of the field induced polarization 
with increasing number of fc-points. Further details will 
be given in Sec. IIVBI 

For the structural characterization in the absence of 
an external electric field, the atoms were allowed to relax 
until the maximum component of the force on any atom 
was smaller than 0.01 eV/A, and the maximum compo- 
nent of the stress tensor was smaller than 0.0001 eV/A^. 
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TABLE I. Reference configuration, cutofT radii, and matching radius between the full core charge and the partial core charge 
for the pseudopotentials used in our study. Units in bohr. 



Tc Si O 

Reference 5s^ 5p'*, 5rf°, 4/° 3s^ 3p^ 3d°, 4/° 2s^ 2/, 3d°, 4/ 

Core radius s 2.00 1.77 1.15 

p 2.00 1.96 1.15 

d 3.00 2.11 1.15 

/ 3.00 2.11 1.15 

Matching radius NLCC 1.30 1.50 1.17 

Scalar relativistic? yes no no 



III. STRUCTURAL PROPERTIES OF THE 
CRYSTALLINE PHASES 

1. a-Si02 cristohalite 

a-Si02 cristobalite crystallizes in the same space group 
(P4i2i2, no. 92) and with the same independent 
atomic positions as paratellurite, a-Te02. Two atoms 
are independent by symmetry: one Si atom at position 
(m,u,0), and one O atom at position (x,y,z). The unit 
cell contains four formula units (twelve atoms). In the 
case of a-cristobalitc, the Si04 entities are almost regular 
tetrahedra (see Fig. [1]) and the Si-0 distances are all 
close to 1.60 A (see Table-HJ). As shown in Fig. [TJ the 
tetrahedra are organized as to form an helical chain along 
the z-direction. While in the directions x and y, the 
structure is made by zig-zag chains. 



(a) 



(b) 



•V 





■ .J,, 



FIG. 1. (Color online) Schematic view of the Q:-Si02 cristo- 
balite unit cell in different perspectives. Si and O atoms are 
represented by grey and red balls, respectively. Solid lines 
mark the unit cell, (a) highlights the helical chains formed by 
the tetrahedra along z, while in (b) the zig-zag chains in the 
X direction are clearly observed. 



Theoretical lattice parameters and Wickoff positions 
are reported in Table-HIl together with some experimen- 
tal values for comparison. Although the data summa- 
rized in Tablc-lnl include results obtained with different 
implementations of the density functional theory (differ- 
ences in the electrons included explicitly in the calcula- 



tion, with and without pseudopotentials, different basis 
sets, different ways of sampling the Brillouin zone), a gen- 
eral trend that can be observed is that both PBE-GGA 
and B3LYP hybrid functional yield an overestimation of 
the lattice constant of up to 3 % with respect to the ex- 
perimental values. The Wyckoff positions obtained with 
the different DFT-methods are in very good agreement 
between them [the maximum difference being in the 0{y) 
position], and arc perfectly comparable with the experi- 
mental results. Indeed, this good performance in a tra- 
ditionally complicated system like Si02 (very sensible to 
many approximations, in particular to basis sets^) vali- 
dates the Siesta basis sets and pseudopotentials used in 
the present work. 



2. a-Te02 phase 

Q;-Te02 crystallizes in a tetragonal unit cell with the 
space group symmetry P4i2i2 {Df, no. 92) as discussed 
beforei^ An schematic view of the unit cell is depicted 
in Fig. m 



(a) 



(b) 



Mi 



FIG. 2. (Color online) Schematic view of the a-Te02 unit 
cell in different perspectives. Te and O atoms are represented 
by grey and red balls, respectively, while the blue spheres are 
the lone pairs of the Te atoms. Solid lines mark the unit cell. 
Meaning of the panels as in Fig. [T] 



The tellurium atom occupies the center of triangle 
bipyramids whose basis is formed by two oxygen atoms 
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TABLE II. Lattice constants (in A), and Wyckoff structural parameters for a cristobalite Si02 (space group P4i2i2). PW 
stands for a plane wave method with pseudopotentials. d(Si-0)i and d(Si-0)2 represent the Si-O bond lengths inside the 
tetrahedra. PBE stands for the Perdew-Burke and Ernzerho^ generalized gradient functional, and B3LYP stands for a three- 
parameter hybrid functional (including part of the exact HE exchange.—). For the O-Si-O angle the average between the four 
possible values is shown. 



^ Reference 
^ Reference 
Reference 
Reference 
Reference 



xc-functional 


Siesta 
PBE 


PW^ 

PBE 


All electron ^ All electron ^ 
B3LYP B3LYP 


Expt. ^ 


Expt. ^ 






Cell parameters (A] 








a 


4.994 


5.073 


4.989 


5.028 


4.983 


4.957 


c 


6.936 


7.085 


6.902 


7.013 


6.955 


6.890 








Atomic positions 








Dl(U ) 


U.oUo 


0.300 


0.307 


0.299 


U.oUu 


U.oUo 


0{x) 


0.236 


0.238 


0.236 


0.240 


0.251 


0.238 


0{y) 


0.118 


0.108 


0.119 


0.104 


0.095 


0.111 


0{z) 


0.186 


0.182 


0.186 


0.178 


0.156 


0.183 








Bond lengths (A) 








d(Si-0)i 


1.629 


1.646 


1.629 


1.615 


1.535 


1.600 


d(Si-0)2 


1.638 


1.646 


1.632 


1.626 


1.606 


1.620 








Angles (deg) 








Si-O-Si 


142.00 


144.39 


142.28 


146.58 


158.86 


144.55 


< 0-Si-O > 


109.32 


109.73 


109.78 


109.75 


119.96 


109.39 



and by the tellurium electronic lone pair, and whose 
apexes are also oxygen atoms. Therefore, the Te atoms 
is coordinated with four O atoms. This Te04 bypi- 
ramidal unit, building block of the tellurium oxides dis- 
cussed in the present work, is referred to as a disphenod. 
Two different Tc-0 bond can be distinguished within 
the disphenod, being the equatorial O atoms closer to 
Tc than the apical O atoms (experimental distances of 
1.87 A and 2.12 A, respectively). As in a-Si02 cristo- 
balite, the polyhedra are connected by vertices to form 
a three dimensional network, ressembling to an helical 
chain along z direction and zig-zag chains along x and y 
directions (sec Fig. [2j) 

Unit cell lattice parameters and internal coordinates 
are reported in Table-HIIl The structural parameters ob- 
tained with Siesta are in very good agreement with those 
obtained with a plane wave code with ultrasoft pseudopo- 
tentials and an energy cutoff of 30 Ry, showing the good 
performance of the basis set used in the present work. 
As usual, the standard overstimation of the experimen- 
tal equilibrium volume by the generalized gradient ap- 
proximation is found. The calculated independent bond- 
lengths at the theoretical equilibrium lattice parameters 
also overstimates the experimental numbers. Neverthe- 
less, the difference between the short equatorial and the 
long axial Tc-0 bond lengths is preserved within Siesta. 



3. 'y-TeOi phase 

7-Te02 crystallizes in an orthorhombic unit cell with 
the space group P2i2i2i (DI, no. 19)..^ A schematic 
view of the unit cell is represented in Fig. |31 This phase 
is metastable at normal conditions, and has been re- 
cently identified by x-ray powder diffraction of recrystal- 
lized amorphous Te02 doped with oxides. The unit cell 
contains four formula units (twelve atoms), with three 
atoms independent by symmetry: one Te atom located at 
(u,u,ti;), and two oxygen atoms labeled as 0/ and O//. As 
in the a-phase, the structure of the 7-phase can be con- 
sidered as a 3D network of corner-sharing Te04 disphen- 
ods. However, in the 7-phase the disphenods are strongly 
deformed, so the length of the four Te-0 bonds are rather 
different (experimentally the bond lengths range from 
1.86 A to 2.20 A), with a much larger spread (0.34 A) 
than in the a-phasc (0.24 A). Indeed, if we assume that 
the longest Te-0 distance (marked with a dashed line in 
Fig. [3]) is too long to form a chemical covalent bond, then 
the structure can be viewed as an infinite zig-zag chain 
of TeOs units in the z direction, connected by the bridge 
Tc-0//-Te. Including now the longest bond Te-0, the 
disphenoids Te04 forms zig-zag chains along x direction 
and helical chains along y direction but with different 
bridges (Te-0/-Te and Te-0//-Te) (see Fig. H) 
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TABLE III. Lattice constants and Wyckoff structural pa- 
rameters for paratellurite a-Te02 (space group P4i2i2, Df, 
no. 92). PW stands for a plane wave calculation performed 
with the Quantum-Espresso packagei^ Oeq and Oap repre- 
sent, respectively, the equatorial and apical oxygens within 
the disphenods. Both Siesta and plane wave simulations 
have been carried out within the PBE-GGA functional.— 
All electron results were computed with the B3LYP hybrid 
functional]^ Units of the lattice constants and distances in A 
and angles in degrees. 

Siesta PW^ AU electron^ Expt. ^ 

Cell parameters (A) 

a 4.987 4.990 4.899 4.808 

c 7.606 7.546 7.792 7.612 

Atomic positions 

0.0261 0.0272 0.0276 0.0268 

0.1418 0.1467 0.1389 0.1368 

0.2494 0.2482 0.2585 0.2576 

0.1973 0.1968 0.1845 0.1862 

Bond lengths (A) 

1.955 1.944 1.909 1.879 

2.192 2.118 2.160 2.121 

Angles (deg) 

104.6 103.6 103.2 103.4 

169.8 171.2 168.1 168.0 

136.2 137.1 139.1 138.6 



Te(?i) 
0{x) 
0{y) 
0(z) 

d(Te-Oeq) 

d(Te-Oap) 



Ocq Tg Ocq 



Te-O-Te 



Reference 
Reference 
Reference 



In Fig. [3l despite the fact that the bridge between Te 
atoms along the z-orientcd chains is always through a 
O// atom, the length of the Tc-0// bond is different, 
ranging between 1.94 A for one of the bonds of the chain 
to 2.02 A for the second bond. 

The theoretical lattice parameters and independent po- 
sitions of the atoms are reported in Table- ITVl The good 
agreement between Siesta and plane waves results con- 
firms and highlights the transferability of our basis set, 
that was optimized for the a-Te02 structure. Again, 
the PBE-GGA functional overstimates the experimental 
equilibrium volume, although the deviation in this case 
(9 %) is slightly larger than usual. This overstimation 
translates also in a slight overestimation of the spread of 
the bond lengths in the case of Siesta (0.40 A). 

Regarding the first-principles simulations on the struc- 
ture of the Te02 phases, we can summarize that the 
disphenoidal configuration of the Te04 entities are re- 
spected both in the a and 7 phases. The lone pair steri- 
cal effect is thus conserved in our geometry optimization. 
The good comparison between our structural parameters 
and the ones obtained with plane waves^ support the 
use of the numerical atomic orbital method implemented 
in Siesta in the present study. 




FIG. 3. (Color online) Schematic view of the 7-Te02 unit 
cell. Te and O atoms are represented by grey and red balls, 
respectively, while the blue spheres are the lone pairs of the 
Te atoms. Solid lines mark the unit cell, repeated along the 
z-direction for the sake of clarity. The yellow dashed lines are 
the longest Te-O distance. 



(b) 





FIG. 4. (Color online) Schematic view of the 7-Te02 unit 
cell in different perspectives. Te and O atoms are represented 
by grey and red balls, respectively, while the blue spheres are 
the lone pairs of the Te atoms. Solid lines mark the unit 
cell. The two symmetry inequivalent O atoms are labeled as 
Oi and 0/7 . The yellow dashed lines are the longest Te-0 
distance, (a) and (b) highlight the zig-zag chain along the z 
and X direction, respectively, (c) focus on the helical chain 
along the y direction. 



IV. NONLINEAR OPTICAL PROPERTIES 
A. Methodology 

When an intense light (a powerful laser) goes through 
an insulator, this medium responds nonlinearly. Its po- 
larization can be expressed as a Taylor expansion of 
the electric field £. 
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TABLE IV. Lattice constants and Wyckoflt structural pa- 
rameters for 7-Te02 (space group P2i2i2i, D2, no. 19). 
PW stands for a plane wave calculation performed with the 
Quantum-Espresso packagei^ d(Te-0)i and d(Te-0)2 stand 
for the short and the long equatorial Te-O distance, while 
d(Te-O)-^' and d(Te-0)2' represent the long and the short ax- 
ial bond lengths within the disphenod. Both Siesta and plane 
wave simulations have been carried out within the PBE-GGA 
functionali^ Units of the lattice constants and distances in 

A. 



where Eq is the permittivity of free space. Unless oth- 
erwise is specified, we will use the SI system of units 
throughout the paper. In the same way, the coefficients 
of the quadratic and cubic terms are used to define the 
second-order and third-order nonlinear susceptibilities, 



(2) 



1 d^P, 
2eo dS.dSk 



(3) 



Siesta PW^ Expt. 



Cell parameters (A) 



a 


5.181 


5.176 


4.898 


h 


8.636 


8.797 


8.576 


c 


4.446 


4.467 


4.351 


Atomic positions 




Te(?i) 


0.9632 


0.9581 


0.9686 


Te{v) 


0.1034 


0.1032 


0.1016 


Te{w) 


0.1368 


0.1184 


0.1358 


Oi(x) 


0.7770 


0.7641 


0.759 


Oiiy) 


0.2935 


0.2851 


0.281 


Oi{z) 


0.1750 


0.1645 


0.173 


Oii{x) 


0.8616 


0.8599 


0.855 


Oiiiy) 


0.0380 


0.0406 


0.036 


07/(2) 


0.7259 


0.7131 


0.727 


Bond lengths (A) 




d(Te-0)i 


1.912 


1.900 


1.859 


d(Te-0)2 


1.983 


1.960 


1.949 


d(Te-0)i/ 


2.116 


2.119 


2.019 


d(Te-0)2' 


2.315 


2.252 


2.198 



Reference 
Reference 



3 



2 dEjdSk 



^ 1 d^V, 



j,k,l=l 



6 dSjdSkdSi 



(1) 



where fc and I refer to cartesian directions, and Vf 
is the zero-field (spontaneous) polarization. The coeffi- 
cients of the previous expansion, derivatives of the po- 
larization with respect to the electric field of increasing 
order, are related with the electrical susceptibilities of 
the material. In particular, the coefficient of the linear 
term is directly proportional to the linear dielectric sus- 
ceptibility (a second-order rank tensor). 



(1) 

At J 



1 dV, 

So d£j ' 



(2) 



^(3) 



1 



6eo dSjdSkSi 



(4) 



Replacing Eqs. ([2])- (01) h^to Eq. ([T|), then the expansion 
of the polarization can be written as 



3 = 1 j.k=l 
3 

+ ^oxf^jli^j^kA + ■ ■ ■ 
j,k,i=i 



(5) 



In general, the polarization depends both on the va- 
lence electrons and the ions. In the present work, we deal 
only with the electronic contribution. The main reason 
behind this approximation is that the second (SHG) and 
third (THG) harmonic generation experiments leading to 
the second and third order susceptibilities are performed 
at frequencies high-enough to get rid of ionic relaxations 
but low enough to avoid electronic excitations 42i This 
constraint implies that both the frequency of £ and its 
second and third harmonics are lower than the fundamen- 
tal adsorption gap. Indeed, SHG and THG experiments 
are pump-probe settings which arc sensitive only to the 
electronic contributions. 

For a completely anisotropic crystal, Eq. ([5]) implies 
respectively 9, 27 and 81 elements for x^^\ X*'^'' X^^"^ 
tensors. These numbers can be reduced significantly by 
considering the specific crystal class.— The nonvanishing 
elements for the symmetry groups of a-Te02 and 7-Te02 
can be found in Table-|V] and Table TVTl respectively. 

Thus, combining the symmetry- allowed components of 
the susceptibility tensors with a judicious choice of the 
direction of the applied electric field, we can restrict the 
expansion in Eq. ([5]) so that only a given component of 
the susceptibility is present. Then, its value can be ob- 
tained by fitting the polarization versus electric field de- 
pendence as obtained from first-principles computations 
of the response of the bulk materials to an homogeneous 
electric fieldJ^iS 

To better fix this procedure, let us develop the way the 

(2) 

Xyxz of 7-Te02 is computed. Taking into account the 
nonvanishing optical susceptibilities of its crystal class, 
P2i2i2i, the expansion of the Eq. ([5]) in this case is 
reduced to 



7 



TABLE V. Symmetry-allowed components of the linear di- 
electric susceptibility, x'^'i ^-nd the second-order, x'^'i 
third-order, x'^'j nonlinear optical susceptibilities tensors for 
the P4i2i2 space group, in which both the Q-Te02 and the 
Q!-Si02 cristoballite crystallize. 







XX 


xyz — — 


yxz 


yy 


xzy — — 


yzx 


zz 


zxy — — 


zyx 



.(3) 



xxxx — yyyy 

zzzz 
yyzz — xxzz 
yzzy — xzzx 
xxyy = yyxx 
yzyz — xzxz 
xyxy = yxyx 
zzyy — zzxx 
zyyz — zxxz 
zyzy — zxzx 
xyyx = yxxy 



TABLE VL Same as in Table-|V]but for the P2i2i2i space 
group, in which the 7-Te02 phase crystallizes. 







^(3) 


XX 


xyz,yxz 


xxxx, yyyy, zzzz, yyzz, xxzz, yzzy, xzzx 


yy 


xzy, yzx 


xxyy, yyxx, yzyz, xzxz, xyxy, yxyx, zzyy 


zz 


zxy, zyx 


zzxx, zyyz, zxxz, zyzy, zxzx, xyyx, yxxy 



- y —i-y r coXyy <--y 



+ f y(2) £ £ ^ y(2) £ £ 

A- Fn (\-'-^'> £ £ £ + v(^) £ £ £ + v(^) £ £ £ 

^ yXyyyy'^y'^y'^y I Ayyzz'^y'^z<^z ^ Xyzzy'-'z'-'z'^y 

A~ Fn f v(^' £ £ £ A- v(^) £ £ £ + v(^) £ £ £ 

^ "--O yXyyxx'-'y'-'x'-'x ~ Xyzyz'-'z'^y'^z ^ Ayxyx'-'x'-'y'-'i 

+ eo (x^y'IxySxS.Sy) ■ (6) 



(2) 

To isolate the Xyxz component, we can apply a field with 
only X and z components, £ = {£, 0, £), so the expansion 
in Eq. ([6]) reduces to 



-Py = n 



(7) 



Assuming the Klcinman symmetry conditions;^ that 
states that the nonlinear optical susceptibility tensor, as 
defined in Eq. ([3]), is symmetric under a permutation of 

the k indices so Xyxz = Xyzx, then Eq. ([7]) transforms 
into 



Up to now we have apphed only basic definitions 
and symmetry properties. The ingredient from first- 
principles comes from the computation of the field- 
induced electronic polarization (the ionic cores are 
clamped at the equilibrium zero-field positions) as a func- 
tion of the external electric field. Here, we have used the 
method of Refs. [isl and [l^ In these milestone works, 
Souza, Kiguez and Vanderbilt showed how to compute 
stationary states of a periodic system at a finite, con- 
stant electric field £ on a uniform discrete fc-point mesh, 
providing that the magnitude of the field does not ex- 
ceed a critical value £c{N), that decreases as the num- 
ber of fc-points N increases. The algorithm, described in 
Sec. V of Ref. [l^ is based on the diagonalization of a 
field-dependent Hermitian operator, Tkl^) in the nota- 
tion of Ref. |T§, for every /c-point in the mesh. The field- 
dependent operator depends explicitly on a given k, but 
implicitly couples neighboring fc-points. Therefore, to 
know the occupied Bloch-like cigenstates, the diagonal- 
ization has to be iterated until the procedure converges 
at all the fc-points and the occupied subspace stabilizes. 
Once self-consistency has been achieved, the overlap ma- 
trix between the cell periodic parts of the Bloch functions 
at neighboring fc-points can be used to compute the com- 
ponent of the polarization parallel to a vector of the re- 
ciprocal lattice, as usual in the discretized formulation of 
the polarization of a solid as a Berry-phaseJ^ 

Once a set of polarization versus electric field data have 
been obtained from first principles, one can make a choice 
between different methods to extract a value for the non- 
linear susceptibility. The first one is a direct fit of the 
raw data (see the symbols of Fig. [5]) to expressions like 
Eq. ©. 

The second one is to compute the derivatives of the 
macroscopic polarization with respect the electric field 
by finite differences, using the Richardson's extrapola- 
tion to estimate the limit ft, — > from calculations with 
two different step sizes: 

/(")(x) = ^D(-\x,h) - il?(«)(x,2ft) + 0(ft4), (9) 

where D^^^ is given by the centered finite difference ex- 
presion 



2h 



£)(2) 13 given by 



(10) 



(11) 

and D'-^' is given by 
(3) ^,_ f{x + 2h) - 2 fix + h) + 2 fix -h)- fix - 2h) 



Vy = + 2eox|,'U'- 



(8) 



=/ ix)+Oih^). 



(12) 
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2,0x10" 




FIG. 5. Field-induced polarization along y, Vy, as a func- 
tion of an applied electric field only with components along 
the X and z directions, £ = {£,0,£), for the 7-Te02 phase. 
The x-ax.is represents the component of the field along the x- 
direction. Atomic coordinates and unit cell lattice parameters 
are clamped to the optimized structure at zero external field, 
and only the electronic response is considered while comput- 
ing the polarization. Circles are the first-principles results 
and the solid line is the fit to Eq. ([8]). 



Here, in the results quoted below, we have used a field 
step ofh = 0.04 V/A. 

Comparison between the results obtained with the two 
methods will be presented in the next Sec. IIVBI 



B. Convergence studies 

It is well known that the total energy ground state cal- 
culations for insulators converge expeoncntially fast with 
respect to fc-point sampling. However, while a given k- 
point sample might be perfectly adequate for some prop- 
erties, it might constitute too severe approximation for 
others 4^ In other words, the convergence with respect 
the fc-point sampling is property-dependent. That is the 
case of the computation of the polarization and its deriva- 
tives, when a discretized Berry phase polarization ex- 
pression is used 4^ In the formalism developed in Ref. 
[l9l . and summarized in Sec. IIV Al both the calculation 
of the stationary states of the periodic system at a fi- 
nite constant electric field, and the Berry-phase polar- 
ization is made on a uniform discrete fc-point mesh in 
the first-Brillouin zone, and the convergence for finite 
field simulations is considerably slower than in total en- 
ergy or charge density calculations!^ The situation does 
not improve significantly when a perturbation approach 
is used instead of the finite field method. Previous first- 
principles simulations on the nonlinear optical suscepti- 
bilities in the framework of the density functional pertur- 
bation theory^ show that the convergence of the results 



is quite poor with the number of spetial fc points, at least 
when the discretization of the formula for the Berry phase 
of the occupied bands is performed after the perturbation 
expansion (although the situation improves dramatically 
when the perturbation expansion is performed after the 
discretization). We can wonder how rapidly our results 
converge with respect to the fc-point sampling. 

In Fig. [6] we represent the Xyxz component of the 
second-order NLO susceptibility tensor of 7-Te02 [panel 
(a)], and the x^xxx component of the third-order NLO 
susceptibility tensor of a-Te02 [panel (b)] as a function 
of the size of the TV x iV x iV shifted Monkhorst-Pack^ 
grid. The results are shown for the two different methods 
used to obtain the derivatives of the macroscopic polar- 
ization with respect the electric field: the direct fit to 
the polarization versus field curve and the Richardson 
extrapolation to compute finite differences derivatives. 
The results provided by both methods are in good agree- 
ment, with differences smaller than 2 % for reasonable 
sizes of the Monkhorst-Pack mesh. For the rest of the 
paper all the reported results have been obtained with 
the Richardson extrapolation method. 

In any case, as we can see from Fig. [51 the convergence 
with respect the number of fc-points included in the sim- 
ulations is rather slow. A way of predicting converged 
values^ for a given magnitude p would be to explapolate 
to — > oo through a least square fit against an analyt- 
ical formula of the form p = poa + a/N^ (where Poo, i, 
and b are adjustable parameters). However, this would 
require several calculations at different A'^. Here, we have 
proceed in a different way, computing only the values of 
the NLO susceptibilities for A^ ~ 20. This would lead to 
results with errors of usually around 2% for the second- 
order and up to 15 % for the third-order susceptibilities. 
However, a similar trend is expected for estimations of 
the susceptibilities in different phases with approximately 
the same unit cell volume and number of atoms per unit 
cell. Therefore, for the consequence of the main goal of 
this work, that is, the comparison of the nonlinear optical 
susceptibilities between phases to ascertain the origin of 
their large values, the previous approach is justified. 

V. RESULTS 

A. Second order NLO susceptibilities 

The values for the symmetry allowed components of 
the second order susceptibilities of a-Te02 and 7-Te02 
are gathered in Table- TVlIl where they are presented in 
terms of the d tensor, often used in the literature of non- 
linear optics, and defined as 

da = d,jk = ^X?jl, (13) 

where the first index, i, refers to a cartesian direction 
(1 for X, 2 for y, and 3 for z), while the second index I 
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jk 



FIG. 6. (a) Xyai^z component of the second-order NLO sus- 
ceptibility tensor of 7-Te02, and (b) the x^^xx component 
of the third-order NLO susceptibihty tensor of a-Te02 as a 
function of the size of the N x N x N shifted Monkhorst- 
Pack^ grid. Circles and solid lines represent the results ob- 
tained with a direct fit of the polarization versus electric field 
curve. Squares and dashed lines have been obtained with the 
Richardson extrapolation formula to compute derivatives by 
finite differences. The lines are fits to analytical functions 
that would allow an extrapolation to A'^ — > cx). 



contracts the two other cartesian indices j and k in Voigt 
notation (see Table IVlII) ) . 



Material Element Theory Experiment 



a-Te02 


di4 


0.0 






0.0 




dse 


0.0 


7-Te02 


di4 


lf.6 




^26 


11.6 




^36 


11.5 



1.4 1.9 



Reference 4^ 
^ Reference li^ 

TABLE VII. Symmetry allowed values for the second order 
nonlinear susceptibilities of a-Te02 and 7-Te02. Units in 
10"^ esu (1 esu = 4.192 xlO"* m/V). 



According to (i) the conditions imposed by the space 
groups (see Tables |V] and IVl)) . and (ii) the Kleinman 

(2) 

symmetry, which means that xlik symmetric under a 
permutation of i, j, and k, all the elements of the NLO 
susceptibility tensor of a-Si02 and a-Te02 should vanish, 
and all the independent elements of 7-Te02 should be 
equal. To illustrate this, let us take, for instance, the du 
element for q:-Tc02. The crystal symmetry imposes that 



1 XX 

2 yy 

3 zz 

4 yz = zy 

5 zx — xz 

6 xy = yx 

TABLE VIII. Relationship between contracted index I and 
cartesian indices jk in the definition of the d tensor in Eq. 
113. 



Il4 



iv(2) 



2 Aya:; 



^25 J 



while following the Kleinman symmetry 



1^(2) ^ iy(2) 

2 ^xyz 2^ /^yxz 



^25- 



(14) 



(15) 



The only way to satisfy Eqs. ([Ti|) and (|15l) simultane- 
ously is fii4 = ^25 = 0. A similar reasoning for 7-TCO2, 
where in principle all the crystal symmetry allowed ele- 
ments of the second order nonlinear susceptibility might 
be independent, shows that the Kleinman rule forces 
them to be equal. A permutation of indices impose that 



dxyz dxzy dyxz dy^x dzxy 



J-zyx- 



(16) 



Both results are clearly visible in Table IVIII 

Unfortunately, the comparison with the experiment 
is not straightforward. The only experimental results 
published up to now on crystals concerns the a-Si02 
critoballite,— and Q;-Te02 phase,— resulting in a weak 
SHG response (indicating the presence of a second order 
susceptibility) with the SHG efficiency of Q;-Te02 five 
times larger than in Si02. The metastable nature of the 
7-Te02 phase makes impossible to grow single crystals 
big enough to be optically characterized. 

It can be surprising at first sight to obtain x'^^ val- 
ues for the a-Te02 where, as explained before, the com- 
bination of the space group symmetry and the Klein- 
man's rules^ should inhibit the presence of a SHG sig- 
nal. Different mechanisms have been called to explain 
this apparent inconsistency. First, Kleinman's relations 
are expected to breakdown if the second harmonic fre- 
quency is close to an absorption band of the material, 
but this was not the case for Q;-Te02 under the mea- 
surement conditions!^ Second, although the Kleinman's 
rule is always satisfied in non-dispersive media, it can 
be broken in dispersive materials, as might be the case 
for the experimental Q;-Te02 samples. As the present 
calculations were done in the static finite field limit, no 
dispersive effect can be taken into account so that the 
Kleinman's conditions are verified. 
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B. Third order NLO susceptibilities 

A glance to Tables |V] and IVII reveals that there are 
21 elements of the third order nonlinear susceptibility 
allowed by spatial symmetry considerations. A first- 
principles estimation of all of the independent elements 
would be a rather cumbersome task, that is out of the 
scope of this work. Instead, we have computed the xxxx, 
yyyy, and zzzz elements of the third order suscepti- 
bilities, as they are informative enough for the struc- 
ture/properties relationship considerations. For those 
particular components both the polarization and the elec- 
tric fields involved are directed along the same cartesian 
direction [see Eq. (jj])]. Results for Q;-Si02, Q;-Te02 and 
7-Te02 are summarized in Table- HXl 

The results obtained with Siesta using the variational 
implementation of the finite field in periodic system ap- 
proach are in very good agreement with those obtained 
with Crystal06, where the field was introduced as a 
sawtooth potential. However, the last approach requires 
the use of a supercell to adapt the periodicity of the po- 
tential and, therefore, the increase of the computational 
cost of the simulation. 

Unfortunately, the comparison with the experiment is 
not so straightforward, since there are no measurements 
yet concerning the crystalline phases of Te02 or Si02. 
Therefore, the values reported in Table- ITXl can be consid- 
ered as purely predictive. The only x^'^' available data for 
tellurium oxides has been measured on the correspond- 
ing glass.— However, Raman spectroscopy measurements^ 
have shown that the structure of 7-TCO2 is close to the 
structure of this glass, so it is tempting to think that 
the order of magnitude of the third-order susceptibility 
should be the same in both the glass and the crystalline 
7-Te02 phase. 

As can be drawn from Tablc- HXl the third order sus- 
ceptibilities of Te02 and Si02 crystals obtained from our 
first-principles simulations arc of the same order of mag- 
nitude as the experimental values for the relevant glasses. 
The very large x'^-* predicted for the two tellurium oxide 
polymorphs, on the range of 10~^^ esu, are nearly two or- 
der of magnitudes larger than for crystallized silica, with 
an average x^^^ (Te02) / x^'^^ (Si02) ratio close to 50, 
thus extending the experimental observations on glasses 
to the case of crystalline compounds. 

We can wonder now about the origin of the large val- 
ues for x'-'^-' in Te02. Despite the fact that Q!-Si02 and 
Q;-Te02 crystallizes in the same space group, their third 
order susceptibilities arc different by two orders of mag- 
nitude. This fact points out that structural arrangement 
by itself can not be responsible for the remarkable NLO 
properties of Te02. 

A more detailed analysis of the x^'^'' tensor elements 
reveals that the highest value are for the z direction in 
the case of a-Si02 and a-Te02, and for the y direction 
for the 7-Te02. As stated in Sec. IHIi those are pre- 
cisely the crystalline directions where the chains display 
an helical shape. Moreover, for the two structures that 



Material Element of x'^' 



This work All electron^ Expt. 



Q-Si02 

Q-Te02 
7-Te02 



xxxx 
zzzz 
xxxx 
zzzz 
xxxx 

yyyy 

zzzz 



0.4 
0.7 
17.3 
31.3 
12.2 
23.6 
20.8 



0.52 
0.59 
18.36 
32.07 



0.28 



14.1 



^ Reference 
^ Reference 



TABLE IX. Symmetry allowed values for the third order 
nonlinear susceptibilities of a-Si02 cristoballite, a-Te02 and 
7-Te02. The all electron simulations have been carried out 
with the B3LYP functional^ as implemented in the Crys- 
tal06 package.— The experimental values, written in italics, 
correspond to the corresponding glass phase, whose structure 
for the Te oxide is not so different from the 7-Te02 phase as 
suggested by Raman spectroscopy measurementsi^ Units in 



10" 



esu (1 esu = 1.398 xlO" 



7v^ 



are structurally similar (a-Si02 and Q!-Te02) nearly the 
same ratio x^zzzz/x^xxxx of around 1.8 is obtained. Since 
no lone pair is present in Si02 compounds, Te lone pair 
effect (orientation of the LP with respect to the electric 
field for example) cannot be responsible for these large 
variations. 

This suggests that high x^'^'' values are structurally in- 
duced and that helical chains are much more favourable 
than the zig-zag chains structures shown along the 
X direction for these materials. Indeed, Mirgorodsky 
et alA and Soulis and coworkersi^ have shown how 
there is a strong link between the structure of poly- 
mer TCO2 molecules and its nonlinear optical properties, 
with the chain-like species developing a drastic increase 
in their third-order hypcrpolarizability with increasing 
chain length. The Te-O-Tc bridges play a dominant role 
in the polarization properties of the long Te02 chains. 
This was attributed to an exceptionally strong nonlocal- 
ity of the electronic polarization, that is, assuming that 
the electric field applied at a given point would induce 
a dipole moment not only at the very point but in the 
vicinity of this point (extending up to eight-ten neighbors 
from the point of perturbation). 

Now, let us turn our atention to the zig-zag chains of 
the two Te02 compounds (along the x and y direction 
for the q;-Tc02 phase and along the x and z direction for 
the 7-Te02 phase.) They are all very similar in shape 
and made of asymmetric single Te-O-Te bridges with the 
bond length sequence —S — L — S ~ L— (where S and L 
stand for short and long bond lengths, respectively). We 
can define an asymmetry index as AI = 100 x {L — S)/L, 
whose value amounts to 11 for the chains parallel to x or 
y in a-Te02, 15 parallel to x in 7-Te02, and 4 parallel 
to z in 7-Te02. It is interesting to note that the AI 
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values evolve as the inverse of the x'-'^^ values, suggesting 
that, for a given chain, the more symmetrical the bridge 
the higher the third order nonUnear optical susceptibility 
value. 



VI. CONCLUSION 

The second and third order nonlinear optical suscepti- 
bilities of two bulk crystalline phases of Te02, and a-Si02 
cristoballitc have been computed using the variational 
approach to compute the response of a periodic solid to 
an external electric field. 

The third order nonlinear susceptibilities are in good 
agreement with previous more expensive theoretical pre- 
dictions, were the electric field is introduced by means of 
a sawtooth potential, and with the experimental results 
for related glass phases. Indeed we were could reproduce 
the large values for the x^^'' tensor elements of the tel- 
lurium oxides, 50 times larger than in pure silica glasses. 

Our results provide some clues to explain the origin 
for the high hypersusceptibilities and the large variations 
with respect to the crystalline directions. In particular 
these properties could be attributed to the presence of 



helical chains in the structure. 

A next step to be taken in order to through some light 
into the origin of the large values of the nonlinear optical 
susceptibilities would require the calculation of the center 
of the localized Wannier functions, and their variation 
with the external fields. 

Our results demonstrate how first-principles calcula- 
tions are an efficient tool to estimate nonlinear optical 
susceptibilities of crystalline solids. These might con- 
tribute to fill the gap usually left by experimental mea- 
surements due to the difficulty in growing single crystals 
big enough to be optically characterized. 
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